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Abstract. We present two complementary methods, each applicable in a different 
range, to evaluate the distribution of the lowest eigenvalue of random matrices in a 
Jacobi ensemble. The first method solves an associated Painleve VI nonlinear differ- 

£SJ ' ential equation numerically, with suitable initial conditions that we determine. The 

second method proceeds via constructing the power-series expansion of the Painleve VI 
*\ '■ function. Our results are applied in a forthcoming paper in which we model the dis- 

rf\ ' tribution of the first zero above the central point of elliptic curve //-function families 

of finite conductor and of conjecturally orthogonal symmetry. 
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^ ■ 1. Introduction 

b , 

We present techniques for calculating numerically the distribution of the lowest eigen- 
value (or synonymously, we say the 'first eigenvalue') of random matrices in Jacobi en- 
sembles JJy' . We proceed as follows. We introduce the Jacobi ensemble Jjy = J-jy' 
o£ N x N random matrices. We relate the distribution of the lowest eigenvalue of ma- 
trices in Jjv to the probability E^' (0; I) that a Jacobi ensemble has no levels in some 
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interval I = [t, 1] for < t < 1. We use two complementary methods to evaluate E^ ' 
relying on its interpretation as the Okamoto r-function of a Painleve VI system, along 
with an auxiliary Hamiltonian function h{t) for which Forrester and Witte [5] have es- 
tablished explicit differential equations of Painleve VI type. Our first method uses the 
Selberg-Aomoto integral to obtain explicit initial conditions for the Painleve IV equa- 
tion satisfied by h(t) which are valid close to the edge t = 1. With these in hand we 
provide the MATLAB code (relying on its built-in ordinary differential equation solver) 

to numerically evaluate h(t) alongside E^ . The second method uses the Painleve VI 
equation with explicit initial conditions at the edge t = together with power series 
manipulations to recursively find the power series expansions of h(t) and E^ . We im- 
plement this algorithm on SAGE using its ability to perform power series manipulations 
and symbolic algebra. 

The use of these two complementary methods is essential in order to compute h(t) 

and E^' (t) accurately over the whole range < t < 1. The Painleve VI equation and 
its solutions have singularities at the edges t = and t = 1. The first method uses a 
solution found starting from an explicit initial condition at a point to = 1 — e close to 1, 
where e > is a small positive parameter we determine empirically. Such an explicit 
initial condition is found in Sections 0] and [5J however, the initial condition is correct 
only up to terms of size 0(e 2 ). The errors introduced by such approximation and by the 
numerical Runge-Kutta method result in a computed solution whose range of reliability 
may not extend to t close to the singularity at t = 0. The second method, described in 
Section [UJ constructs a truncated but otherwise exact power series for h(t) about t = 
(the singularity at t = is handled indirectly) up to terms of order 0{t D+1 ) where D 
is the degree of the truncation. Such a solution is reliable over any interval [0, u] with 
u < 1, provided D is large enough, though not necessarily over the entire interval [0, 1] 
in view of the singularity at t = 1. In Section [7] we analyze the range of parameters 
a, b, N for which both methods are stable in the sense that both numerically computed 
solutions agree in some subinterval [u, v] of (0, 1), which implies that the numerical solver 
is robust for this range of parameters. 

It is important to note that the methods we use apply to non-integer values of N. 

Painleve differential equations have played a role in many problems in random matrix 
theory, ranging from the distribution of the eigenvalues in the bulk to the largest and 
smallest eigenvalues, and have been extensively studied; for our purposes, the most 
relevant are the investigations of solutions to Painleve VI. We briefly mention some of 
the literature. We refer the reader to the special edition of the Journal of Physics A 
(Volume 39, Number 39, 2006), which celebrates 100 years of Painleve VI, especially 
the historical introduction and survey [H] and the article by Forrester and Witte [6] on 
connections with random matrix theory; see also the recent works by Dai and Zhang |2j 
and Chen and Zhang pQ for determinantal formulas obtained from ladder operators. 

The main contribution of this paper is the derivation of an algorithm to compute 
numerically the distribution of the lowest eigenvalue in the Jacobi ensembles, and a 
discussion of its implementation and accuracy. The motivation for this project comes 
from attempts to understand the observations in [T3] on the distribution of the first zero 
above the central point in families of elliptic curve L-functions when the conductors 
are small. The Katz-Sarnak conjectures 0(2] predict that as the conductors of the 
elliptic curves tend to infinity their zero statistics should agree with the N — > cxd scaling 
limits of the corresponding statistics of the eigenvalues of matrices from a classical 
compact group. For suitable test functions this was proved in |12|I17]: however, for finite 
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conductors the numerical data in |13j is in sharp disagreement with the limiting behavior 
of these random matrix ensembles. In particular, the first zero above the central point 
is repelled, with the repulsion decreasing as the conductors increase. In a forthcoming 
paper we complete the study of the low lying zeros of elliptic curve L- functions, and 
obtain a model which describes the behavior of these zeros for finite conductors. One 
of the key ingredients in our model is the lowest eigenvalue of these Jacobi ensembles 
of N x N matrices, often requiring non-integer values of N, which is the main result of 
this paper. 

2. Jacobi ensembles and their first eigenvalue 

Let Jjv = Jjy' denote the Jacobi ensemble on N levels < Xj < 1, j = 1,2, . . . ,N, 
with real parameters a,b > — 1. Explicitly, the TV-level (joint) probability density of 
levels of Jn on [0, 1]^ with respect to its Lebesgue measure dx\ dx2 ■ ■ ■ dxjy is given by 

N 

(2.1) Cp b) HW( Xj ) J] {x k - X] f 

where the weight function W = W^ a,b ' on [0, 1] is given by 

(2.2) W(x) = x b (l - x) a 

and Cjy' is the ensemble's normalization constant. Jacobi ensembles as described 
above correspond to suitable ensembles of self-dual random matrices via the angular 
variables cpj defined by 

1 + cos ( 
Y 

Note that the edges x = 0, x = 1 correspond respectively to (p = ir, cj) = under this 
change of variables. We refer the reader to [3] and the forthcoming book [7] for details 
regarding the matrix realizations of Jacobi ensembles, for which we will otherwise have 
no direct use. In what follows we will go back and forth between the abscissas Xj and 
the angular variables (f)j, but will in any case refer to the associated ensemble by the 
Jacobi name and denote it by Jn- In terms of the angular variables, and with respect 
to Lebesgue measure on (0, n) N , the TV-level (joint) probability density for J^ is given 
by 

JV 

(2.4) C { x' P) JJiy(cos^) \\ (cos0 fc -cos^) 2 

j=l l<j<k<N 

with the weight function w = w^ a '^' on (0, n), 

(2.5) w(cos (p) = (1 - cos <j)) a {l + cos <pf '. 

The parameters a, /3 > — \ are related to a, b above by a = a + 1/2, /3 = b + 1/2, 

and Cjy-' is the appropriate normalization constant, namely CL' = 2~'" +a+fe+2 )" 

C^ for the constant C^ of (|2.ip . For suitable choices for a and @ we obtain 

the joint probability density of the iV independent eigenphases for the classical groups 
of matrices SO(2iV), SO(2iV + 1) and USp(2iV), when the latter are endowed with an 
invariant (Haar) probability measure and regarded as random matrix ensembles. The 
case a = f3 = corresponds to SO(2iV), a = 1 and (3 = corresponds to SO(2iV + 1), 
and a = ft = 1 to USp(2iV). This is explained in detail in [3]. Below in (|4.23[) we give an 



(2.3) Xj = 7r — 1 , < 4>j < it. 
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explicit expression for the normalization constant Cjy ' • ( Jacobi-distributed pseudoran- 
dom sequences of levels can be generated from a uniform pseudorandom sequence using 
only the Jacobi joint probability density via for instance the Accept-Reject Algorithm 
whose applicability is quite broad; see for instance [To].) 

As remarked above, the Jacobi ensemble Jn describes the eigenvalue statistics in 
suitable ensembles of self-dual random matrices having N pairs of eigenvalues e ±4 ^, 
j = 1, 2, . . . , N; we call (j> the eigenphase of the eigenvalue e 1 ^ . Let E^ ' (n; I) denote 
the probability that a random matrix A G J/v has exactly n eigenphases in the inter- 
val / = [0,0]. As shorthand notation we write E^' (4>) f° r ^n (0; [0, </>]), namely 
the probability of having no eigenphases in the interval [0,0]. The probability density 
function u^' ((f)) of the distribution of the first eigenphase is related to E^' ((f)) by 

(2-6) 4^((f)) = -±E^\(f)). 

We can deduce the relation (|2.6p as follows: assume that the interval [0, (f)} contains no 
eigenvalues. Then a small increment e > of the interval to [0, + e] has two possible 
outcomes. Either the interval [0, (f) + e\ contains no eigenvalues, or it contains some. The 
probability of the first event is E^' P) \(f) + e). It follows that E^' P) ((f)) - E { ^ ] \<j> + e) 
is the probability that the interval [(f), <fi + e] contains at least one eigenvalue; as e — > 
there can be only one, namely the first eigenphase in [(f>, (f> + e\. Thus 

(2.7) lim grWz£TW±f) _ -'&>«) 

e-^0 £ d(f> 

indeed yields the probability density function v£ (<f>) of the first eigenphase. An 

alternative way to prove (|2.6[) is to observe that 1 — E^' ((f)) is the probability that 
[0, (p] contains at least one eigenphase, hence that the first eigenphase 4> m i n is at most (f>; 
otherwise said 1 — Ej^' '((f)) is the cumulative distribution function of the first eigenphase 

0miru so its derivative is equal to the probability density function v£ (</>)• In Section [5] 
we shall need to scale the angular variable by a factor of N/tt in order to consider 
eigenphases of mean unit spacing on [0, N]. 
We have 

e { ^\<p) = cp p) f r - f f[(i - cos^ra + cos^y 

(2.8) H H i =1 

x j J (cos (f>j — cos (f>k) d(f)\ ■ ■ ■ dipN 

l<j<k<N 

for fixed a, f3 > —1/2 and the normalization constant C^' of (|2.4p . There is no known 

method to evaluate the multiple integral in equation (|2.8p exactly. E^' ((f)) is related 
to a Painleve VI transcendental function h(t), namely a certain solution to a second- 
order nonlinear ordinary differential equation. In Proposition 14.41 we provide the first 
few terms of a power-series expansion of E\r' ((f)) for (j) close to 0; these provide the 
initial conditions for the differential equation we aim to solve. Our main reference for 
the theory is the work of Forrester and Witte [5] . Their result is stated for the abscissal 
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counterpart to the function E*? '((f)) of (|2.8|) . namely the function E^' ' (i) defined by 

/•* /■< N 

(2.9) E$ h >(t)=c$ b) •••/ n^ 1 -^) n (^-^) 2 ^i---^iv 

JO JO - =1 l<j<k<N 

with the normalization constant Cjy' of (|2.ip . The functions (|2,9p and (|2.8p are related 
by the change of variables 

(2.10) a = a - 1/2 and b = (3 - 1/2 
along with 

(2.11) t = l+|os0 ^ ^1 + cos^ 

where < t, Xj < 1. Explicitly, 

(2.12) E^\ 4) ) = E^( l + C °^ 



2 

3. First method: The auxiliary Hamiltonian and Painleve VI 

Both of our mutually complementary methods rely on the interpretation of E^' 
as an Okamoto r-function and ensuing relation to a Painleve system with associated 
auxiliary Hamiltonian h(t); this Hamiltonian arises as the solution of a Painleve VI 
equation with the exact parameters determined by Forrester and Witte in Proposition 
13 of [5] as follows. 

Proposition 3.1. Let a,b > —1 and N be a positive integer. The auxiliary Hamiltonian 

(3.1) h(t) = t ■ e' 2 [b] - ±e 2 [b] + t(t - 1)| log^' 6) (t) 

where 

,, , , , % /, T a + b a — b a + b , T a + b 

(3.2) b = (61,62,63,64)= 1^ + — >— ' 2 _ '" 7V "^~ 

(3.3) e 2 [b] = 6i6 2 + 6163 + 6164 + 6263 + 6264 + 6364, 

(3.4) e' 2 [b]= 6163 + 6164 + 6364, 

satisfies the following Painleve VI equation in Jimbo-Miwa- Okamoto a -form 

4 

(3.5) h'it) (t(l - t)h"{t)) 2 + (h'(t)[2h - (2t - l)ti(t)} + hb 2 b 3 b 4 ) 2 = ]J(ti(t) + b 2 k ). 

fc=l 

Furthermore, we have the boundary condition (as t — > 0) 

(3.6) m = (-i e2[b] - *(* + «)) + (4m + w(Ar+ ^ + ° +t) ) ' + °« 2 >- 

Note that besides simplifying the notation in Proposition 13 of [5] we also swap the 
a's and 6's therein. The parameter a is equal to the order of vanishing of the Jacobi level 
density at the edge t = 1, whereas a + 1/2 is the order of vanishing of eigenphase density 
at the edge <fi = 0. We remark that the apostrophe in the symbol e 2 has no specific 
meaning and is merely used to visually distinguish it from e 2 (in a manner consistent 
with the notation of reference [5]), whereas the apostrophe in h! and elsewhere in this 
manuscript means differentiation: ti(t) = *§, h"(t) = ^, E^' b) '(t) = f t E^' b) (t), etc. 



DUENEZ, HUYNH, KEATING, MILLER, AND SNAITH 



Pay close attention to the fact that the initial condition given in (|3.6p holds at t = 0. 
This condition will be used in Section [6] to construct a power-series solution. For our 
intended application, however, we are most interested in the behaviour of E^ (</>) for <p 
close to zero, which in view of the change of variables (|2.1ip corresponds to t close to 1. 
Unfortunately, the singularity of the Painleve equation at t = 1 significantly complicates 
the numerical evaluation of the function h{t) in this range. 

Our first method of solution will numerically compute h{t) starting instead from an 
initial condition given at some fixed point t = to close to 1, say to = 1 — e for some small 
positive e to be chosen empirically. The determination of this suitable initial condition 
is a delicate issue that depends on the analysis carried out in Section |H 

Following Edelman and Persson [3] , we seek to compute simultaneously E^ ' (t) and 
the Hamiltonian h(t) via a (non-autonomous) differential equation for the triple of func- 
tions 



(3.7) 

of the form 

(3.8) 

with initial conditions 

(3.9) 



H(t) 



EX> u> (t) 



' ? (a,b) 
J N ' 
h(t) 

h'(t) 



Hn 



dH 
~~dt 



H(t ) 



Ft(H) 



EX' u) (t ) 



;(a,6). 
' J N \ 
h(t ) 

h'(to) 



where to = 1 — £ f° r small e > 0. (The singularity of the Painleve equation and its 
solutions at to = 1 preclude taking simply e = 0.) 
From (13.11) we obtain 



(3.10) 



Jt EN {t) 



h{t)-te' 2 [h] + \e 2 [h] 
t(t-l) 



E 



(a,b) 

N 



(t), 



and likewise from (13.5 



(3.11) h"(t) 



1 U j= i(h'(t) + bj) - (h'(t)[2h - (2* - l)h'(t)} + 6i6 2 6 3 6 4 ) 



t(l-t)V 

Therefore, f|3.8[) holds with 
(3.12) 

/ 

fhi(t)\ 



Ft 



h 2 (t) 

\hs(t)J 



h'(t) 

/t 2 (t)-te 2 [b] + ^e 2 [b] 
t{t - 1) 

h z {t) 



hi(t) 



1 ULMt) + 6?) - (h 3 (t)[2h 2 (t) - (2t - i)h 3 (t)} + hb 2 b 3 b A y 



\t(i-t)V h 2 {t) 

The MATLAB code to compute F t (H) is given in the appendix and is also available for 
download at http : //www . maths . bris . ac . uk/~mancs/publicat ion s . html[ We em- 
ploy the built-in ordinary differential solver ode45 from MATLAB, which implements a 



/ 
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Runge-Kutta method giving an approximate solution of (|3.8|) (and thus of the sought 
density z/jy (4>) = — -jxE^ ((f)) of the distribution of the first eigenphase). It remains 
still to determine the initial condition Hq = H(to) for (|3.7|) as follows. We shall find the 
small-e asymptotic behavior of E^' (1 — s) (equivalently, what we will actually do is 
find the small-^> asymptotics of E^' ((f))). By differentiation we then find 

(3.13) pN b) V and W^ h){t) 

and thus (through its definition (|3.ip ) we obtain asymptotically good approximations 
to h(to) and h'(to) for any to close to 1. This gives the triple Hq of initial conditions 
for (|3.7p . In the following section we compute the asymptotic behavior of E^' '((f)) for 
(j) close to 0. 

4. Taylor series expansion for E^f (<f>) 

In this section we compute the probability Ej^' '((f)) that a random matrix from a 
Jacobi ensemble Jn as defined in Section [2] has no eigenphase in the interval [0, (j)] for 
small (f> > 0. As described at the end of Section [31 this enables us to derive the initial 
conditions for the system of differential equations (|3.8p which gives the distribution of 
the first eigenphase (|2.6p . Forrester and Witte [5] consider this same limit in their 
equation (1.38), but here we derive a further term in the approximation. Our result is 
stated in Proposition 14.41 

We require some notation. For n = 1, . . . , N and a,/3> —1/2 we define the integral 

I(n) := tifr® f ••■ f [*... /^(l-cos^ra + cos^)' 3 
Jo Jo Jo Jo „•_., 
(d -n v v 

V^'-V N—n times n times 

X j | (COS (f>j — COS 0fc) d<^l • • • d(f>N. 

Then we have the following lemmata: 
Lemma 4.1. For E^ ((f)) as given in 112. 8\) and for I(n) as defined in f^.ip we have 



J N 

(1.2) <'-"(o| = l-N-lQ.)+r:\I(2) -[")!$) + ■■ ■ + [-!)» I{N). 



*'%) =■ l-N.I(l)+(f)l(2)-r\r,, ,v; 



Lemma 4.2. For 1(1) as defined in |^.1[ ) we /\aue 

_(N-1)H 2 +(§ + §) ft 

w;/zere 

^ ._ r(a + 7V + l/2)r(a + /3 + iV) 






j.2a+3 



and 

(4.5) H 2 



2 2 «r(a + l/2)r(a + 3/2)T(N + l)r(/3 + N - 1/2) 
r(a + N + l/2)r(a + /3 + TV + 1) 



2 2a + 1 r(iV + l)r(/3 + TV - l/2)r(a + l/2)r(a + 5/2) ' 
Lemma 4.3. For /(n) as defined in ft4.1\ ) we have for n > 2 and a > —1/2 

(4.6) I(n) < 02an+2n 2 -n_ 
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We postpone the proof of the above lemmata for a moment in order to state the 
desired Taylor series expansion of i% (4>) for small 4> > 0. 

Proposition 4.4. For E^' (eft) as given in I12.8\) we have 



(4.7) N y ' V 2a + 1 

+ O(0 2Q+4 ) 
where H\ is defined in |^.^|) and 77 2 in H_ 



(^- | )/:/, + (§ + f)i:/i 



k2a+3 



2a + 3 



Proof. Let n > 2. We can apply Lemma 14.31 to every term I(n) in (J4.2J) of Lemma 14. 11 
Thus, I(n) <C (j) 2an + 2n ~ n . As n > 2, we have 2cm + 2n 2 — n > 2a + 4 for every 
q; > —1/2, hence 7(n) <C <^> 2a+4 . Thus every I(n) in Lemma 14. II can be absorbed into 
the error term 0(4> 2a+i ) of 7(1) in ([O]) of Lemma [OJ D 

Now we prove Lemmata 14. 1[ 14.21 and 14.31 

Proof of Lemma \4-l\ Recall that E^ ' ((f)) is the probability of having no eigenphases 
in the interval [0, 4>}. We denote the event that there is at least one eigenphase in [0, <f>] 
by B. Then its complement CS is the event of having no eigenphases in [0, (/>]. Hence 

(4.8) E^ $) (</>) = P(CB) = l-P(B). 
We now focus on P(B). Let 

(4.9) B k : = {(</>!,. ..,</> k ,...,<f> N )€ [0,ir} N where 4> k € [0,0]}, 

which is the event that 4> k lies in [0, <fi] and the remaining eigenphases lie anywhere in 
[0,7r]. Note that Bj and B k ,j / k, are not necessarily disjoint. 

We can write the event of having at least one eigenphase in [0, <p] as 

N 

(4.10) B = |J S fc . 

fc=i 

For the probability of the event Bk we have 

iV 

(4.11) P(B k ) = C { ^> / H(l - cos^) Q (l + cos^) /3 J] (cos^ - cos0 fe ) 

For any fe and any j ^ k we have 

(4.12) P(P fe ) = 7(1) and P(B J n B k ) = 7(2), 
and in general, for 1 < i\ < ■ ■ ■ < i k < iV, 

(4.13) Pftn-nBij = 7(fe). 

By the inclusion-exclusion principle and the symmetry above, 



(4.14) 



:J .i ■ 



p \ 


U B * hLH) 2. p(t^ n---nP 4 




Vfc=l / fc=l l<ii<---<i fc <Af 


N 


P(Pi) - (^) P(Pi n P 2 ) + ^ P(Si n P 2 n P 3 ) 


+ 


■ ■ + (-if+ 1 F(B 1 n-nB JV ). 
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Thus 

(4.15) P(B) = N -I(l)-(^y(2)+(^y(3) + --- + (-l) N+1 I(N). 

Finally, the probability of having no eigenphases in [0, cp] is given by 

N 



P$B) = Cp P) f ■■■ T \\(1 - cos<^) a (l + cos&y 



J. — wo y^7 } V x _r ~ KjKJO V? 
0=1 

x I j (cos 4>j — cos 4>k) dcpi ■ ■ ■ d(/)N 
l<j<k<N 
= 1 - P(B) 

(4.16) = 1 - N- 1(1) + (Vl(2) - (Til® + ■■■ + (~l) N I(N)- 



Proof of Lemma \4.ty First we recall Selberg's integral (see Chapter 17 of |llj): 

af^l-s,)"" 1 

i<j<k<N 



•*° ^° i=l Kj<k< 



_ l 2 7 



u , 7) = ^ r(i + 7 + n)F( P + j 7 )r(>7 + j 7 ) 

1 ' j fi r(i + 7 )r(p + r, + [Ar+i-i] 7 ) 

and Aomoto's extension of Selberg's integral for 1 < R < J\f 

I dx x -.. ( dxxflxjflx^il-x^- 1 H \( Xj -x k )\ 2 ^ 

"'° ^° 3=1 Z=l l<j<k<M 

, 4181 A p + (AA-j) 7 - A p r 1 r(i + 7 + J7)r(p + J7)r(T ? + j 7 ) 

l ' J ^ P + V + (2AA-J- 1) 7 1 = 1 r(l+ j)T(p + V + [M + j - l]j) 

both valid for integer J\f and complex p, rj, 7 with 

(4.19) Re(p) > 0, Re(r?) > 0, Re( 7 ) > -min ( 1- ■ E( ' (/,) H< " ' 7i 



^'(AA-l)'(AA-l) 
The version of Selberg's integral we are interested in is related to (|4.17j) by 
p = s+1/2, f] = r + 1/2, 7 = 1, 

and the change of variables 

1 + cos d>j 
Vj = o 



D 
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as follows (note T(2) = 1): 

/ d(/)i--- d^JJ(l -cos0/) r (l + cos^) s ][ (cos (pj - cos 4> k ) 2 

J° ^° 1=1 l<j<k<Af 

(4.20) = ^ +S - J ) f d yi ... f dy M ([(1 - yi y-^yr 1 ' 2 \{ («j " Vk 

J ° J ° 1=1 Kj<k<Af 



= 2 mN+r + s-i) tt r(2 + j)T(s + 1/2 +j)r(r + 1/2 +j) 
W r(s + r + AA + j) 

The version of Aomoto's extension of our interest has parameters 

7 = 1, 



p = r + l/2, 


rj = s + l/2, 


and we change variables 






1 — cos (hj 

Zj = - 

3 2 


in (|4.18|) to obtain 




(4.21) 




/ #!•••/ d^T\{l- 
Jo JO u-i 


N 

- cos0 fe ) JJ(1 -cos^) r (l 
;— i 



]~J (COS (j)j - COS <j)kY 
l<j<k<N 
1 R N 



Zkf 



= 2 «+^(^+- 1 ) r ^ • • • r cbrv n «* n ^r i/2 (i - -/r i/2 n & 

•^° ^° fe=i z=i i<i<fc<A^ 

= 2 fl +W+r+s -i) A r + l/2 + AA-j ^ r(2 + j)r(r + l/2 + j)r( S + l/2 + j) 
-Mj r + s + 2jV-j [} o T(r + s+J\f + j) 

We can now determine the normalization constant CL' in ()4.1ip . We have 

(4.22) 

AT 

-i /"7T /»7T _ _ 

^N = / d(j)\--- I d4>N TT (1 — COS0j) Q (l + COS (j)jY II (COS (j)j — COS 4>k) 

Jo Jo j=l l<j<k<N 

By setting AT = N,r = a, and s = /3 in ()4.20p we obtain 

(a ov r^r 1 _ N(N+ a +p-i) "tt 1 r(2 + j)r(/3 + 1/2 + j)r( a + 1/2 + j) 
[ } N ~ JLJL r(a + /3 + iV + j) 

Now, for small cp > and using Selberg's integral (|4.20p we wish to evaluate 

(4.24) 1(1) = C£'flj(l) 
where 

1(1) := / •••/ / K(fa,...,<f> N )H(<f> 1 ,...,<l> N )d<t>id(h---d<t> N 

Jo Jo Jo 

with 

JV 

(4.25) K(4> 2 ,...,<t> N ) ■= l[(l - cos ^ (1 + cos cPjf II (cos^-cos0 fe ) 2 , 

3=1 2<j<k<N 
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and 



A? 



(4.26) iJ(0i,...,0iv) := (1-cos(/>i) q (1 + cos0i) /3 ]^[(cos0i -cose 



fc=2 



Now we evaluate H(4>i, . . . ,4>n)- The Taylor expansion around (f>\ = of the first 
factor of H((pi, . . . , <f) N ) in (J4T26D is 



(4.27) 
(1 — cos^i) Q (l + COSc 



A2a ±2a+2 

2 a 2 a ■ 12 



+ 0{4>f 



b 2a 



a 



P 



2 a —fi \ 2 a ~fi ■ 12 2 a— ^+ 2 



2 /3 _ ^~ 2 ^\ + 



2a+2 + o ^2 a+ 4 ) _ 



The Taylor expansion of the terms in the second factor defining H((f>i, . . . , 4>n) m (|4.26p 
is 



n 

fc=2 



COS 01 — COS (pk) 



(4.28) 



T / r/> 2 \ 1 

II 1-f + O(0f) -COS0, 

fc=2 L V - / J 

Af 
[] [(1 - cos0 fc ) 2 - 0f(l - COS^fc) + O(^)] 



fc=2 

Af 

[](l-cos<^) 2 -0i 



fc=2 



Af AT 



3=2 k=2 



(1 - COS^fc) 2 



COS( 



+ o(4) 



Using (|4T27D and (14351) in (ET26D gives 



ff (&,... 



'JV 



(4.29) 



if! 

2 q-/3 

' N 

no 

fc=2 



O 



/? 



2a - /3 . 12 + ^^]^ +2 + o(^ +: '' 



COS^fcJ 



AT JV 

En , 

i=2 fc=2 



(1 -COS^fc) 



COS( 



+ 



Multiplying out and collecting terms by powers of 4>\ gives 



k2a N 



Hi 



>N) 



bW 



2 a ~P 



COS(p k 



J.2Q+2 
,2 01 



(4.30) 



fc=2 



o 



2"-^ 



N N fi ± \2 

y^-pj (1 -COS0fc) Z 

i=2 k=2 ^ 3 



+ 



P 



2 a~f3 . 12 2 a ~/ 3 + 2 



V 



^ +2 n(i-cos^) 2 +o(0 2 



2a+4\ 



fc=2 



12 
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Integration of the last expression (|4.30p gives 



(4.31) 



N 



" 6 2a+1 

if(01,...,0Ar)#i = Y[(l-COS(j) k ) 2 

k=2 \ a + ) 



>a-fi 



' N N C\ J, \2 

ST^TT (1 - cos (j) k y 



3=1 k=1 

N 



COS <pj 



i.2a+3 



(2a + 3)2 



\OL-t 



Hence, to evaluate 1(1) 
to compute 



(4.32) 



(4.33) 



rn pn 



- frn -cos <t> k ) 2 ( — % — + -JL-)f^- + o(<p 2a+5 ). 

k=2 V 7 

= fo " ■ Jo Jo K (<fa> ■■■■> (f>N)H((f>i ,..., 0jv)#i • • • d<f>N, we have 

• d(j)N and 



/ ■■■ / K(<f> 2 ,...,(f>N 

Jo Jo k=2 

N N 



{J (1 -COS</>fc) 2 # 2 



'o jo 



K((f>2,...,(f>N) 



sr^yr (.1 -cos (fa) 

^11 1 - cos (/>,' 

j=2 fc=2 Y -> 



d4>2 ■ ■ ■ d4>N- 



Observe that the integrand of (J4.33P is symmetric in its variables fa , ■ ■ ■ , 4>n ■ There- 
fore we have 



(4.34) 



/o Jo 

= (N -I 



K(fa,... An) 



N N ri j. \2 

(1 - COS0 fc ) Z 

cos 4>j 



En 1 ^ 



h ■ ■ ■ d(p N 



N 



k=3 



P7T p1T ^ 

/ ••• / K((j)2,... ,0at)(1 -COS(f> 2 ) TT(1 - cos (j) k ) 2 dfa ■ ■ 
Jo Jo ,7_q 

Evaluating the integral (|4.32|) yield: 

/ ••• / K((f> 2 ,... ,4>n) TT(1 - cos (f> k ) 2 dfa- --d(j) N 
/o Jo k=2 



■d(j) N . 



(4.35) 



Is 

N 



■ ■ ■ JT (1 — COS (j)j) a+2 (l + COS (f>jf J] (cos <f>j — cos 4> k ) 2 dfa ■ • ■ d(j)N 

•*° j=1 1<j<k<N 

and using Selberg's integral ()4.20p with AT = N — 1, r = a + 2, and s = f3 gives 



cos 

N-2 



>fc) 2 #2 • • • #JV 



(4.36) /*.■. flCfa,..., <f> N )f[(l - 

Jo Jo k=2 

= 2 (N-i)(N + a + fS) TT £(2 + j)T{fi + 1/2 + j)r(q + 5/2 + j) 
=o T{a + P + N + l+j) 
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Normalizing the last expression (|4.36j) with Cjy from (|4.22|) we obtain 

(4.37) Cp P) I ■■■ I K (<f> 2 ,..., <f, N )f[(l- cos fa) 2 : # 2 ---# JV 
Jo Jo k=2 

1 r(a + N + l/2)r(a + /3 + N) 



2 a +P r(a + l/2)r(a + 3/2)F(N + l)T(/3 + iV - 1/2) ' 
We now evaluate the integral in (|4.34j) . For this we note that 

N N N 

(4.38) (l-cos0 2 ) J|(l-cos^) 2 = JJU-cos^na-cos^). 



So 



tZ — O rv — Ai tZ — O 



COS(/>fc) dcf>2 ■ ■ ■ d(f>N 
k=3 

N N 



(4.39) 



Jo Jo k=3 

= •••/ K((f)2,...,(f)N) TT(l-COS0 fc )TT(l-COS0 fc )d0 2 

■ / ° J ° k=2 k=3 

= f ... f JT(1 -008^(1 + cos hf TT ( 
JO JO ~-o o^tt^r 



j=2 2<j<k<N 

N M 

x IJ(1 -cos<^> fc ) Y[(l ~ cos (f) k )d(f>2 ■■■ d(f) N 

k=2 k=3 

r r N N 

= •••/ n( i - cos ^)TT( i - cos ^) a+1 ( i + cos ^) /3 

J ° J ° fc=3 j=2 

X J I (COS <j)j — COS 4>k) d<p2 ■ ■ ■ d(j)N- 
2<j<k<N 

With R = N - 2,J\f = N - l,r = a + 1, and s = /3 in Aomoto's integral (j4T2l"|) this 
evaluates to 

(4.40) 



iV 
I ■■ I K((f>2,...,<t>N)(l-COS<j> 2 )i[(l 

Jo Jo r , 



COS(j)k) 2 d(t)2 ■ ■ -d(j)N 



fc=3 

JV-2 



2 (A r -2)+(Ar-l)(Af+a+^-l) TT « + -/V" + 1/2 - j 

11 a + /3 + 2iV-l-i 
j=i 

^ r(2 + j)r(q + 3/2 + fflty + 1/2 + j) 

2 (N-2) + (N-i)(N +a+ p-i) r(« + JV + l/2) r(a + (3 + N + l) r(a + N- 1/2) 

T(a + 5/2) r(a + /3 + 2iV - 1) T(a + 1/2) 

^ T(2 + j)T(a + 3/2 + j)T(/3 + 1/2 + j) 
1 = 1 T(a + f3 + N + j) 
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Normalizing the last expression (|4,4U[) with C^ from (|4.22|) we obtain 

C N' P) I ■■■ J K(<fo,..., <f> N )(l- COS 2 ) 11(1 - COS ^) 2 #1---#7V 

(4.41) ° ° k=3 

_ 2 - a - -i r(q + N + 1/2)I> + + JV + 1) 

r(7v + i)r(/3 + TV - i/2)r(a + i/2)r(a + 5/2) ' 

Putting (|3~57|) and (JUT]) into (J4T2lj) we obtain 



t2a+l _ A2a+3 



J (!) = ^i ,o ■ 7v^f -(JV-l)g2 



(4.42) 

where 

(4.43) #1 
and 

(4.44) F 2 



(2a + 1)2°-^ v ' (2a + 3)2°-^ 



2 Q -^-12 2 a ~P+ 2 J 2a + 3 

1 r(a + TV + l/2)r(a + $ + N) 

2 a +P T(a + l/2)r(a + 3/2)r(iV + l)r(/3 + N- 1/2) 

1 r(a + iV + l/2)r(a + p + N + 1) 



2 a+f3+i T ( N + 1 )p( /3 + N _ i/ 2 )r(a + l/2)r(a + 5/2) ' 
Finally, we rewrite the result slightly and get 



cb 2a+1 , s (t> 2a+3 fa p\ cj) 2a+3 

1(1) = H 1 ^-—-(N-l)H 2 ^-—-[- + ^)H 1 ^—- + 0(<fy 



(4.45) 



2a + 1 v ' 2a + 3 V 12 4 / 2a + 3 



i2a+5\ 



Hi; 



2a + 1 



(N-l)H 2 +(^ + ^)H 1 



,2a+3 

- + O(0 2q + 5 ) 



where 

(4.46) ffi 
and 

(4.47) H 2 



2a + 3 

r(a + JV + l/2)r(a + p + JV) 

2 2 «r(a + l/2)r(a + 3/2)r(7V + l)r(/3 + TV - 1/2) 

r(a + N + l/2)r(a + p + iV + 1) 



2 2a + 1 r(iV + l)r(/3 + iV - l/2)r(a + l/2)r(a + 5/2) ' 
This completes the proof of Lemma 14.21 



□ 



Proof of Lemma \4-3\ The definition of I(n) is 

Tin) ■- fT" / ••• /" / ••• f TTO -ros,,,)"(.l -c(*^)" 



^ / ••• / / ••• / ri( i - cos ^) a ( i - cos ' 

Jo Jo Jo Jo .•_-, 



(^.4loJ N—n times n times 

X 1 I (COS (j)j — COS (/>fc) d(f>i ■ ■ ■ d(j)N. 

l<j<k<N 

Here we are only interested in the size of I(n) in terms of n, so we can disregard the 
normalization constant Cjy* . For n > 2 we consider 

J(n) := C^ _1 /(n). 
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Then we have 



(4.49) 

I(n) = \ dcj) n+ i--- \ dc/>N IT (I- cos (pj) a (I + cos 4>j 

x J J (cos 4>j — cos4>k) 2 

n+l<j<k<N 

N 
x / < (1 — cos (-.),,)"(. I --cose),,!" I | (cos(fi„ — cos d)j) 



/ < (1 ~COS(fi n ) a (l + COS (fin) 13 |~[ (COS <f> n - COS 
^ ^ j=n+l 

r<t> N 

x / (l-cos^ n _ 1 ) Q (l + cos0 n _ 1 )' 3 ][ 



COS (f>n-l — COS </>■,) 



x (cos(/> n _i — cos0 n ) <i0. 



2, 



'n-1 



AT 

2 



x / (l-cos0i) a (l + cos0i) /3 ff (c 

x (cos 0i — cos 02 ) x (cos (fi\ — cos 03) x • • • x (cos (fil — COS (fin) d(fi\ \d(fi n - 



(COS (pi — COS <pj 
j=n+l 



Now for j, k < n, j ^ k and (fij,(fik G [0, 0] f° r small > we have 



(4.50) (cos0j — cos (fik) < (1 — cos< 



There are (n — 1) + • • • + 1 = n(n — l)/2 terms of the form (cos (fij — cos (fik) (with 
j = 1, . . . , n — 1, and /e = j + 1, . . . , n) occurring in ()4.49p . each of which we bound 
using (|4.50p . 

From equation (|4.3ip in the proof of Lemma 14.21 we derive for k = 1, . . . , n that 



r<t> N 

/ (1 -cos(fik) a (l + cos(j) k )P If (cos (fik -cos (fij) 2 d(j) k 

r , ' /0 3=n+l 

(4-51) 

j.2o+l 

= na-~ ^ 2 24^tt) +0< ^ +3) - 

j=n+l v ^ ; 
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Using (|4.50p and (|4.5ip we obtain 

/ d<f) N 
Jo 

x 



I(n) < i <io ll+ \ 



N Yl (l-cos0 i ) a (l + cos(/ )i ) /3 



j=n+l 
Yl (cos <t>j - cos cf> k f{l - cos <f>) n ( n -V 



n+l<j<k<N 

N 

nc 

j=n+l 



k2a+l 



COS 0, ) X 5- 



A' 



+ 



2a+3\ 



«<^n+l • • • 
^0 



d^iv J[ (1 — cos0j) a (l + cose 



(4.52) 



j=n+l 

x J] (cos^-cos^^^l-cos^™- 1 ) 

n+l<j<fc<iV 

N 



X 



(2a+l)n 



n(l - cosc/> 7 ) 2n x t — ^7 xi- + Off 
v Y3) 2 a ~' 3 (2a + l)} n vr 

=n+l 

AT 



n(2a+l)+2\ 



^n+1 



COSc 



\a+2n , 



1 + cos<^-) /3 



.7=n+l 



X Yl ( COS & - COS <t>kf{l ~ COS 4>) n{ - n ~^ 

n+l<j<k<N 

0(2a+l)n 



[2«-/ 3 (2a + l) 



+ 



n(2a+l)+2\ 



Hence for some constants c\ and C2 



j(n) < ci(l — cos< 



\n(n— 1) 



^(2a+l)n + ^n(2a+l)+2\ 



(4.53) 



c 2 [</» 2 + o{4> 

C-2 



4x-i n(n~l) 



^(2a+l)n + Q^n(2a+l)+2 

(2a+l)n , n (in(2a+l)+2> 



^2n(n-l) + /^2n(n-l)+2^ 

C 2 (/' (2a+1)n+2n(n ~ 1) + 0(^( 2 «+ 1 )™+ 2n ( n - 1 )+ 2 ) 
i2cm+2n 2 -n , q/ j (2a+l)n+2n(n-l)+2\ 



Hence I(n) <C ^ 2an + 2n n 5 which implies /(n) <C 02an+2n n ag re q U j rec [ to finish the 
proof. □ 

Remark 4.5. It is natural to ask whether we can determine further terms of the Taylor 
expansion of E^' {(f)) than provided in Proposition 14.41 By Lemma 14.11 this requires 
us to determine I(n) with n > 2. However, for I(n) with n > 2 we encounter multiple 
integrals which cannot be dealt with using Selberg's or Aomoto's integral. To our 
knowledge the required generalization of Selberg's integral does not exist. 
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5. First method: Initial conditions for the Painleve equation 

Here we use the results from the previous section to actually state the initial conditions 
(|3.9p for our system of differential equations ()3,8p which gives the distribution v^' ((f)) of 

the first eigenphase (|2.6|) of random matrices in the Jacobi ensemble Jn = J^ ■ With 
these initial conditions we then implement a MATLAB algorithm to compute a numerical 
approximation for v^ ((f)). The full MATLAB code is provided in the appendix. For 
its implementation we follow some of the ideas in Edelman and Persson [3] . 

As outlined at the end of Section [3] we are left to provide E^ (to), h(to), and h'(to), 
namely the initial conditions (|3.9|) . for some to = 1 — £ with small e > 0. Recall that, 
according to Proposition 14.41 we have 
(5.1) 



E N 
with 
(5.2) 

and 
(5.3) 



(a,/3), 



1-NlH 



i-i 



h 2a+l 



2a + 1 



12 4 ' 



ffi 



(N - 1)H 2 + 
T(a + N + l/2)r(a + f3 + N) 



k2a+3 



2a + 3 



+ 0(^ a+4 ), 



Ho 



2 2a T(a + l/2)r(a + 3/2)F(N + l)T(f3 + N - 1/2) 
T(a + N + l/2)r(a + (3 + N + 1) 



2 2a + 1 T(N + l)r(/3 + N - l/2)r(a + l/2)T(a + 5/2) ' 

Via the substitution <fi = cos~ 1 (2t — 1), equation (|5.ip provides a good approximation 

for Ejy' (to). In view of the definition (|3.ip of the auxiliary Hamiltonian h, we need 

to differentiate E^' (t) twice with respect to t in order to obtain the initial conditions 
h(to), h'(to) (the second and third components of (|3.9p ). 

Since (b = cos _1 (2t — 1) implies d(f>/dt - 
chain rule and equation (|3.ip give 



1/^(1 -t) and E%> u >(t) 



p(a,b) 



E^ P) ((f)), the 



/i(i) = t • e' 2 [b] - ±e 2 [b] + t(< - ]/ 



d p( a , b )(+\ 
dt h N W 






£!T J (t) 



(5.4) 
and 



£ 



fc'(t) = e' 2 [b] + 



t-e^l-ieatbl + ^l-t) 



l-2t ^' /3)/ (cos- 1 (2t-l)) 
2 v /t(l^i4 Q '/ 3 )( cc . s -i(2t-l)) 



^ / (cos- 1 (2t-l)) 



E { ^\cos- l (2t-l)) 



E^" (cos-\2t - l)i 



(5.5) 




cos~ 1 (2t-l)) 



^/tO^t)Ep 0) (cos-i(2t-l)¥ 

^ Q ' /3), (cos- 1 (2t-l)) 



2^t(r^)^ ) (cos~ 1 (2t-l)) 
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^' /3)/ '(cos- 1 (2t - 1)) ^• /3) '(cos^ 1 (2t - I)) 2 



+ 



^' /3) (cos- 1 (2t - 1)) ^' /3) (cos- 1 (2t - I)) 2 ' 

With (|5,ip . (|5.4p and (|5.5[) we have the initial conditions (|3.9p for (J3.8J) . As mentioned 
earlier, with these initial conditions we can now implement a MATLAB algorithm to 
compute the distribution of the first eigenvalue of random matrices in the Jacobi en- 
semble. We provide the full MATLAB code in the appendix. Also, it can be obtained 
from the authors or from their web pages, such as 



http : //www . maths . bris . ac . uk/~mancs/publications . html 



6. Second method: Symbolic solution using power series 

In this section we describe an algorithm to compute the power series expansion of 
the Painleve function h(t) at t = 0, leading to the numerical computation of Ej^' '((f)) 
for <j) close to it. It is rather unfortunate (for our intended application) that t = 1 is 
a branch-point singularity of h(t), and as a consequence a power-series expansion of h 
about t = 1 is a Puisseaux series (i. e., a series in fractional powers of t), at least if 
the parameters a, b (and, eventually, N) are rational numbers; for arbitrary values of 
the parameters the situation would be even more complicated. Therefore, we content 
ourselves for the time being with finding the power series expansion about t = 0. 

The idea of the algorithm is very simple: The coefficients ho, h\ of the expansion 
h(t) = ho + h\t + /i2i 2 + . . . are given in equation ()3.6p . These are used to bootstrap 
a recursive search for the higher coefficients h 2 , /13, • • • , regarding each unknown % as 
implicitly defined by the earlier coefficients ho, hi, h 2 , • • • , ^fc-i through the Painleve 
equation (|3.5p . Given the complicated nonlinear nature of the Painleve equation it is not 
immediately obvious that this approach will work in practice, but fortunately it does, 
and each successive coefficient h^ is expressible as a rational function of the previous 
ones, hence ultimately as a rational function of the parameters a, b, N. Once many terms 
are computed, the power series for h(t) can be used to evaluate E^ (t) by solving for 
the latter in equation ()3.ip ; then we use the change of variables t — >• <p and differentiation 

to compute the density v^' (<f>) of the distribution of the first eigenvalue, at least for 
</> relatively close to tt. 

We seek to find the coefficients hk as exact rational numbers; for this reason, we will 
work with exact (truncated) power series with rational coefficients. This approach has 
the enormous advantage that the complicated operations needed to evaluate both sides 
of the Painleve equation (|3.5p introduce no numerical errors at all in the evaluation of 
successive higher coefficients. The price paid is a more expensive calculation compared 
to one done using exclusively floating-point arithmetic. Remarks on the choice of the 
number of needed terms to reach adequate numerical precision follow below in Section [71 

In the appendix we list the code for an implementation of this algorithm in SAGE |16j 
(Software for Algebra and Geometry Exploration), a free and open-source computer 
algebra system, although Maxima [10] plays an important role behind the scenes. The 
Python syntax underlying SAGE is clean and the code listing should prove useful both 
SAGE newcomers and those interested in porting it to other computer algebra systems. 
Line numbers from the SAGE code listing will be referenced below as needed. 

We take the parameters a, b and iV to be (fixed) rational numbers, and regard the 
function 

(6.1) h{t) = h + hit + h 2 t 2 + • • • 
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as having coefficients which are rational numbers. In particular, from equation f|3.6|) we 
have 

(6.2) ^ = _l e2 [ b ]-iV(6 + iV) 

(6-3) h 1 -e 2 [b] + ^—^ 

with e 2 [b], e' 2 [b] given in terms of a,b, N by (pT2"|) . (pT3l) and (|Q]) . 

In the SAGE code listing, the values of the basic parameters a, b, N are hard-coded, 
as well as the maximum degree DEGREE of precision of all (truncated) power series (lines 
1-4). Note that the algorithm's implementation depends crucially on inputting rational 
values for the parameters a, b, N. With a limitation explained below, the algorithm can 
handle integral as well as rational values of the parameter n = N. 

When, at any given point, the coefficients ho,. . . ,hk-i are known, but hk is still 
unknown, we wish to regard the latter as an indeterminate, say hk = X. Let 

(6.4) h(t) =ho + ht+--- + hk^- 1 + Xt k . 

Then h lies in the ring Q[X, t] of polynomials in X, t with rational coefficients. Let us 
denote by LHS(h) and RHS(h) the polynomials obtained by substitution of (|6.4p in 
the Painleve equation (pTKjl . and let PEZ(h) = RHS{h) - LHS(h) ( "Painleve-Equal- 
to-Zero" — lines 33-37). The natural hope is that if ho, h\, . . . , h^-i are chosen correctly, 
then PEZ(h) = pk(X)t + 0(t +1 ) for a non-constant polynomial pk £ Q[a?] an d some 
(unspecified) polynomial 0(t +1 ) divisible by t + . Then hk should be chosen to be a 
root of pk(X). 

Performing exact polynomial arithmetic in Q[X, t] to compute PEZ(K) is quite ex- 
pensive, especially since we only need to determine the coefficient p(X) of the lowest 
power of t. For computational purposes, however, it is enough to regard h as a finite trun- 
cation of an infinite power series in t, systematically neglecting any higher-order terms 
not needed for the immediate purpose at hand — namely the determination of Pk(X). 
Fortunately, SAGE can do algebra in power series rings (with help from Maxima). 

Henceforth we work in the ring S = Q[AT][[£]] of power series in t whose coefficients 
are in Q[X] (polynomials in X with rational coefficients) as done in lines 13-14 of the 
SAGE code. In order to determine Pk(X) it is enough to know h up to 0(t k+2 ) terms. 
(It is not quite enough to work modulo t k+1 because Pk(X) a priori depends also %+i, 
not just on the currently unknown X = hk', luckily, the solution found a posteriori shows 
this dependence to be fictitious.) 

We thus set 

(6.5) h(t) = h + hit + ■ ■ ■ + h k ^t k - 1 + Xt k + 0{t k+2 ) 

and compute PEZ(h) = pk{X)t k + 0(t k+1 ) for a linear polynomial p^ (the only excep- 
tion is pi, which is quadratic with a trivial root X = 0). In the SAGE implementation, 
we store the coefficients hk in a SAGE list g (lines 26 & 27 — it is here that crucial use 
is made of the initial conditions (|6.2p and (|6.3p ) . 

Successively (main loop in lines 31-46), for each k = 2, 3, . . . it suffices to take hk as 
the unique (nontrivial) root of pk(X), which is a rational number. (Note that this is 
the only point in the algorithm at which symbolic algebra is needed, and that because 
of the trivial nature of the equation solved it would be easy to write an implementation 
dispensing with any use of symbolic algebra, a task which we presently avoid simply 
because the resulting code would be longer and more difficult to read.) This root is found 
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in lines 38 & 39, then appended to the coefficient list in line 40 and h(t) reconstructed 
using the newly-found h\~ in lines 42-45. 

A couple of remarks on the code are in order. SAGE does not seem to understand that 
we wish to interpret the variable x of the power series ring F to be a "symbolic" variable 
with respect to which the equation PEZ[i] == (i. e., p%{X) = 0) is to be solved. We 
therefore need explicitly to replace the ring's variable x with a symbolic variable X (lines 
29 & 38) before finding the root of PEZ [i] , which is then appended at the end of the 
coefficient list g (lines 39 & 40). We also found it simpler to reconstruct h from scratch 
in lines 42-45 using the coefficient list g than figuring out a way both to (i) increase 
the precision of h, and (ii) substitute the newly-found rational coefficient g[i] for x. 
(The wrapper QQ( . . . ) around the argument of g is used to convert ("coerce" in SAGE 
lingo) the root found by Maxima to a bona fide SAGE rational number.) 

Solving for E^ (t) in the definition (|3.ip of the auxiliary Hamiltonian h(t) yields 



/ / ■ h(t) - t ■ e' 2 \b] + ±e 2 [b] 
i(i-l) 



E]£' (t) = exp [ j — — — — dt + c ] for a suitable constant c 




where C = e c and /i (1) (t) = (h(i) - h )/t = h\ + h 2 t + h 3 t 2 + ... 




*(* — 1) J t-1 

(6.6) i 

— e 2 [b] - N(N + b) from flX 

- Cexp L(N + i) log, + / ^(t)-4M-N(N + b) \ 

- «"<»+» exp (j mt)-m-w+v> A 

Note that the integrand in the last expression above is regular at t = 0. We define 

(6.7) ^expQ'^-^-^ + ^A 

Then J-(t) has a series expansion (in integral powers) about with J~(0) = 1, and 
equation (|6.6p reads 



(6.8) E ( fi' b \t) = Ct N ( N+ VT(t). 

It follows that the leading-order term of the power series expansion of E\S' (t) about 
t = is Ct . ' ■ ' , explaining the nomenclature leadexp ("leading exponent") for the 
auxiliary SAGE variable defined in line 25. 

The value of the constant C in (|6.8p can be read off from equation (3.5) in [5] as 
the quotient of the normalization constants C^' = 1/Sn(q + 1, b + 1; 1) for the Jacobi 
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ensembles Jj^' and C N ' = 1/«Sjv(1, b + 1; 1) for Jfo' (recall that we swap the role of 
a and b relative to Forrester and Witte; moreover, CJy' = C N ' a ). Explicitly, 

( . q) C^ &) ^(1,6+1;!) 

1 j df 6) 5 w (a+l,6+l;l)- 

The SAGE function Jac(a,b,n) in lines 16-22 evaluates the Selberg integral Sn(cl-\- 
1,6 + 1; 1) for integer values of N = n using formula (|4.17|) : the value of C is then 
computed and stored in leadcoef (line 24). 

This particular implementation naturally depends on JV = n being a positive integer. 
However, it is possible to evaluate Jac(a,b,n) in closed form using Barnes' G-function. 
Unfortunately, the G-function is not yet implemented in SAGE; however, given any 
(future) implementation thereof, the following code 



# *NOT* 


VALID SAGE CODE UNTIL BARNES' 


G 


IS 


IMPLEMENTED 


! ! ! 


def Jac 


(a,b,n): 












return 


G(n+a+b + 2.)/G(2*n+a+b + l.) 
*G(n+a + l.)/G(a + 2.) \ 
*G(n+b + l.)/G(b + 2.) \ 
*G(n + 2.) 


\ 











would correctly compute leadcoef and the program would be capable of handling gen- 
eral rational values of n. Alternatively, the value of leadcoef can be computed by any 
other means and manually input into the SAGE code as a hard constant. 
We now rewrite (|3.10|) in the form 

(MO) <><«> = (^k_£kW _ W+S) -£» w , 

which, together with ()6.8p allows computing E^' (t) and its derivative E^ (t) as done 
in lines 58-64 of the code. Finally, the cumulative distribution function 

l-E^\® = l-E^( 1 -+^ 

and its derivative (cf. (|2.6p and (|6.10p ) 

uW>w = —Efcnw = ±Ie^> f 1 + COS( 



2 JV V 2 



can be computed directly, as done in lines 66-72. (Note that h sincj) = y/t(l — t) if t = 
— 7p^.) The entire SAGE code can be obtained from the authors and is also available 



for download at http://www.maths.bris.ac.uk/~mancs/publications.html 



7. Comparison of the two methods 

As might be expected, our numerical implementation in MATLAB of the Painleve 
solver does not work equally well in all parameter regimes. In this section we describe 
the tests we have carried out on the code and the conclusions about the parameter 
regimes where a robust solution can be obtained. The MATLAB (Runge-Kutta) solver 
starts from initial conditions near <j> = (that is, t = 1) and numerically extends the 
computed solution towards (j> = tt (or 6 = —cj) = N in scaled units). The power series, 
on the other hand, is an expansion around eft = tt (or, equivalently, t = 0); it is therefore 
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accurate at the opposite end of the interval on which we are solving. Hence, if the tail 
of our numerical solution matches the initial behaviour of the series solution, we are 
confident that the numerical solver has worked correctly. 

There are three parameters to vary in the input to the numerical solver: a = a — 1/2, 
b = j3 — 1/2 and N. There are also three variables we can adjust in the MATLAB code 
to try to coax a solution: to, the starting point near t = 1 (<f> = 0) for the numerical 
solver; reltol and abstol, which control the accuracy of the numerical solution. 

After testing the code for various values of iV (integer and non-integer) from about 1 
up to 100, it appears that a good solution can be found on a standard desktop machine in 
a few seconds with to = 1 — 10~ 7 , reltol= 10~ 5 and abstol= 10~ 6 for any —0.5 < a < 
and —0.5 < b < 0.5 (the range for b that is relevant to the classical groups). In Figure 
[1] we see examples of code that runs efficiently and matches the series expansion in the 
tail of the distribution. 



a=0 b-0 N-2 



a— 0.5 b=0.5 N-5 



series solution (100 terms) 
numerical Painleve solver 
initial conditions 



series solution (99 terms) 
numerical Painleve solver 
initial conditions 



Figure 1. Plot of v 



a+l/2,6+1/2 

N ' 



jtO), the scaled distribution of the first 



eigenvalue. On the left N = 2, a = and 6 = 0. The numerical solver was 
given the initial conditions (blue dotted line) and produced the solution 
shown with the dot-dashed line. This is indistinguishable from the series 
expansion (green solid line) using 100 terms. On the right N = 5, a = 
—0.5 and b = 0.5. The numerical solver was given the initial conditions 
(blue dotted line) and produced the solution shown with the dot-dashed 
line. The tail agrees with the series expansion (green solid line) using 
99 terms. In both figures the numerical solver was run with values of 
t = 1 - 10~ 7 , reltol = 10~ 5 and abstol = 10~ 6 . 



For values of a > the MATLAB solver breaks down. Trials with a = 0.001 still 
work, but already a = 0.01 fails to produce a good solution. Moving to away from 1, 
for example to 1 — 10~ 2 , helps achieve a better curve, but as can be seen from Figure 2, 
the initial conditions are not close enough to the true curve at this point to produce a 
valid solution. Decreasing reltol and abstol, even by a factor of 1000, does not make 
a visible difference to the curve. We note that unfortunately this means that while the 
MATLAB solver works very well for the group SO(2A r ), we cannot use it to produce 
solutions for SO{2N + 1) and USp(2iV). 

For a > and small values of N a solution valid over the whole interval can still be 
glued together by matching the series solution (with a sufficient number of coefficients) 
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2.3 
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a=0.5 b— 0.5 N=5 



series solution (300 terms) 

numerical Painleve solver 

initial conditions 



0+1/2,6+1/2 

N 



(jtO), the scaled distribution of the first 



Figure 2. Plot of u 

eigenvalue. On the left N = 2, a = 0.5 and b = —0.5. The numerical 
solver was given the initial conditions (blue dotted line) and produced 
the solution shown with the dot-dashed line. This fails to produce an 
accurate solution, as shown by comparison with the series expansion 
(green solid line) using 50 terms. On the right N = 5, a = 0.5 and 
b = —0.5. The numerical solver was given the initial conditions (blue 
dotted line) and produced the solution shown with the dot-dashed line. 
This fails to produce an accurate solution, as shown by comparison with 
the series expansion (green solid line) using 300 terms. The value of to 
in the two plots are to = 1 — 7 x 10~ 2 and to = 1 — 10 -2 respectively, and 
for both plots reltol = 10 and abstol = 10 . 



for the tail and bulk of the curve with the known asymptotic behavior near = 0. The 
right-hand plot in Figure [2] shows that this is certainly possible for N = 5. 

The series solution produces a very accurate curve with only 50 terms when N = 2, 
but the number of terms needed to obtain a solution which is meaningful over a large 
interval increases with N, which is to be expected because the most interesting behavior 
occurs near cj> = and can only be captured at the price of using many terms in an 
expansion about <p = tt. A good curve for N = 5 requires around 300 terms. We did not 
produce results for higher TV" as the run time was prohibitive, but on a fast computer 
more terms could be computed and good solutions for larger N could be achieved by 
this method. 



8. Summary 

In summary, we find that our MATLAB code for numerically solving the nonlin- 
ear second-order differential equation (Painleve VI) of the auxiliary Hamiltonian h(t) 

associated to the r- function E^' (t), giving the distribution of the first level in a Ja- 

cobi ensemble J 
in the ranges: 

interval [—0.5, 0.5] as this is the range of interest interpolating between the classical 
compact groups SO(2iV), SO(2iV + 1) and USp(2./V). The program's numerical accu- 
racy is confirmed by comparing it with a power-series expansion of the solution found 



N , appears to work fine for arbitrary iV provided the parameters lie 
—0.5 < a < and —0.5 < b < 0.5. We restricted our tests to the 
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using SAGE. Unfortunately the restriction on a to be non-positive means that the nu- 
merical solver cannot cope with the symplectic group USp(2iV) nor the odd orthogonal 
group SO(2iV + 1). For positive a and small values of N this limitation can be over- 
come by using the series expansion to obtain the tail and the bulk of the distribution, 
matching the series result with the initial conditions for the behaviour near the origin. 
In a forthcoming paper our algorithms are applied to study the distribution of the first 
zero of L-functions with even functional equation associated with quadratic twists of a 
fixed elliptic curve. In the limit of large conductors in families of such L-functions the 
zero statistics are expected to be modelled by eigenvalues from SO(2iV). 

Appendix A. MATLAB Code 

Here we give the MATLAB code for numerically solving the Painleve VI equation 
associated to the distribution of the first eigenphase of random matrix ensembles in the 
Jacobi ensemble Jn = Jjy' ■ 

The main program is painleve6.m and it computes the numerical solution to the 

system of differential equations (|3.8j) for the Jacobi ensemble jjy ' . For fixed variables 
a, b, N the program is called by entering 

(A.l) function [t ,H, theta, Fp] = painleve6(a,b,N) 

Here t, H correspond to the variables t and H (see (|3.7p ). and theta corresponds to the 
rescaled angular variable 9 = N(j)/ir = — cos _1 (2t — 1). The output variable Fp in (jA.lj) 
is a vector of values of the rescaled distribution 

(\ o\ (0+1/2,6+1/2) {nd\ it . (■k9\E^ (t) f x , 

(A-2) ^ \^vj=^v sm ^j7(^iy^)-^[b] + ,e 2 [b]] 

of the first eigenphase (the rescaling achieves mean unit spacing of the N eigenphases 
#i, . . . , 0jv on [0, N]), as obtained from (|3,1U|) via the rescaled variable 6 (= theta). The 
subsequent command 

(A. 3) plot (theta, Fp) 

plots the distribution Fp (as defined in (|A.2[) ) of the first rescaled eigenphase 6 m \ n for Jjy. 
Notice that a=-0.5 and b=-0.5 corresponds to SO(2A). Likewise setting a=0.5 and 
b=-0 . 5 gives SO(2iV + 1) and finally USp(2iV) would correspond to choosing a=0 . 5 and 
b=0.5. 

The code of painleve6 .m is given as follows 
function [t ,H,theta,Fp] = painleve6(a,b,N) 

tO = 1 - le-7; 

phiO = acos(2*t0-l) ; 

% in parameter regions where the numerical solver is robust tO can be 
% taken to be about l-le-7 

bl = (a+b)/2+N; 

b2 = (a-b)/2; 

b3 = -(a+b)/2; 

b4 = -(a+b)/2-N; 

e2p = bl*b3 + bl*b4 + b3*b4; 
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e2 = e2p + b2*(bl+b3+b4) ; 

r = a+0.5; % note that r and s correspond to alpha and beta 

s = b+0.5; 

% initial conditions below are from Section 5. The files Hone.m 
% and Htwo.m calculate often-used ratios of gamma functions. EN, 
% ENderivl and ENderiv2.m calculate E_N~{(a,b)}(phi) , and its 
% first and second derivatives with respect to phi 
% (*not* with respect to t) . 

rootO = sqrt(tO*(l-tO)) ; 
EO = EN(r,s,N,phiO); 
EOd = ENderivl (r, s, N, phiO) ; 
EOdd = ENderiv2(r,s,N,phiO) ; 
HO = [ EO; ... 

t0*e2p - 0.5*e2 + rootO*EOd/EO; ... 

e2p + (0.5-t0)/root0 * EOd/EO - EOdd/EO + (E0d/E0)~2 ]; 

% First component of H will be \tilde{E}_N"{(a,b)}(t) 

% Second component of H will be the auxiliary hamiltonion h(t) 

% Third component of H will be h' (t) 

7, HO contains the initial conditions at t=tO (Note: tO is a 

7o number very close to 1, not close to zero! ! 

opts=odeset ( ' reltol ' , le-5 , ' abstol ' , le-6) ; 

% a command like "opts=odeset( 'reltol' , le-5, 'abstol' , le-6) ; " works 

7o fine for parameter ranges where the numerical solver is robust, 

7o and the program just takes a few seconds/minutes to run. 

7o As a becomes positive the differential equation 

7o solver has trouble and we don't get a correct solution 

7o when N is large, a plot on a better scale is produced by 

7o replacing the second argument in the range [t0,0.01] with 

7. 0.5*cos(2*pi/N)+l/2 

[t,H] = ode45(@p6diff , [tO.O.Ol] ,H0,opts,bl,b2,b3,b4,e2p,e2) ; 

F = H(:,l); 

h = H(:,2); 

7o theta is the scaled angular variable: theta=(N/pi)*acos(2*t-l) 

theta = (N/pi)*acos(2*t-l) ; 

7o Fp(theta) is the distribution of the first eigenvalue in scaled 

7o variables 

Fp = (pi/2/N)*sin(pi*theta/N) . *F./t ./(t-1) . *(h-e2p*t+e2/2) ; 

The system ()3.12|) is defined in the file p6dif f .m. The code is given as 
function dH = p6diff (t ,H,bl,b2,b3,b4,e2p,e2) 
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dH = zeros(3, 1) ; 

dH(l) = H(l)/t/(t-l)*(H(2)-e2p*t+e2/2); 
dH(2) = H(3); 
dH(3) = +sqrt( . . . 
( ... 
(H(3)+bl~2)*(H(3)+b2"2)*(H(3)+b3~2)*(H(3)+b4"2) . . . 
- ( H(3)*(2*H(2)-(2*t-l)*H(3)) + bl*b2*b3*b4 )~2 ... 
) / H(3) ... 
) /t/(l-t); 

The often used ratios of gamma functions Hi and H<i (see ()5.2p and (|5.3[) ) are coded 
in the files Hone .m and Htwo .m given as 

function HI = Hone(r,s,N) 

HI = gamma(r + N +0.5) * gamma(r + s + M) / 2~(2*r) ... 

/ gamma(N + 1) / gamma(r + 1.5) / gamma(r + 0.5) ... 

/ gamma(s + N - 0.5) ; 

and 

function H2 = Htwo(r,s,N) 

H2 = gamma(r + N + 0.5) * gamma(r + s + M + 1) / 2"(2*r+ 1) ... 

/ gamma(N + 1) / gamma(r + 2.5) / gamma(r + 0.5) ... 

/ gamma(s + N - 0.5) ; 

Finally, we provide the code for EN.m, EMderivl.m, and ENderiv2.m 

function E = EN(r,s,N,phi) 

% This is an expansion in phi around phi=0 of E_N~{(a,b)}(phi) . 

7, See (ISTTD . 

exponent = 2*r+l; 

E = 1; 

E = E - N*Hone(r,s,N)*phi . "exponent /exponent ; 

exponent = exponent + 2; 

E = E + N*( ... 

(N-l)*Htwo(r ,s,N)*phi . "exponent . . . 
+ (r/12.0+s/4)*Hone(r,s,N)*phi. "exponent ... 
) /exponent ; 



function Ed = ENderivl(r,s,N,phi) 

7, This is \frac{dHd\phi}E_N~{(a,b)}(phi). 

7, See (IB~4D . 

exponent = 2*r; 

Ed = -N*Hone(r,s,N)*phi. "exponent ; 

exponent = exponent + 2; 

Ed = Ed + N* ( ... 

(N-l)*Htwo(r,s,N)*phi. "exponent . . . 
+ (r/12 . 0+s/4) *Hone (r , s , N) *phi . "exponent 
) ; 
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function Edd = ENderiv2(r,s,N,phi) 

I This is \frac{dHd\phi}E_N~{(a,b)}(phi). 

% See (15751) . 

exponent = 2*r-l; 

Edd = -2*r*N*Hone(r,s,N)*phi. "exponent ; 

exponent = exponent+2; 

Edd = Edd ... 

+ N*(exponent+l)*( ... 

(N-l)*Htwo(r,s,N)*phi . "exponent . . . 
+ (r/12.0+s/4)*Hone(r,s,N)*phi. "exponent ... 
); 

Appendix B. SAGE Code 

Below follows the SAGE code implementing the symbolic power series evaluation of 
the r-function. Note that the parameters a,b,N correspond to a, b, N, which along 
with DEGREE (the degree of the sought-after truncation of the power series) are hard- 
coded in the first four lines. Note also that the code provided only works for rational 
values of a, b and integer iV (cf., Section [6j). Once run, the code defines functions E(t) , 
Ep(t) , pcummul(phi) andnu(phi) implementing E^' (t), iE^' (t), 1—E^((p) and 
v(4>), respectively. In particular, nu has been used to produce the plots in figures 1 and 2. 

Note that in Python syntax any text following the literal # is simply a comment. 

i DEGREE = 50 # degree (plus 1) of the truncated power series 

2 a = -1/2 

3 b = -1/2 

4 n = 4 

5 

6 bl = (a+b)/2+n 

7 b2 = (a-b)/2 

s b3 = -(a+b)/2 

9 b4 = -(a+b)/2-n 
10 e2p = bl*b3 + bl*b4 + b3*b4 
n e2 = e2p + b2*(bl+b3+b4) 

12 

13 R.<x> = QQ['x>] # R = Q[x] 

14 S.<t> = PowerSeriesRing(R, default_prec=DEGREE) 

#S=Q[xJ[[tJJ 

16 def Jac (a , b , n ) : 

17 jac=l. 

is for i in xrange(n): 

19 jac *= real (gamma (a+i +1.)) * real (gamma(b+i +1.))\ 

20 * real (gamma( i +2.)) / real ( gamma ( a+b+n+ i +1.)) 

21 return jac 

22 # NOTE: Jac (as above) only works for * integer* values of n 

23 
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24 leadcoef = Jac (0 ,b , n)/ Jac (a ,b , n) 

25 leadexp = n*(n+b) 

26 g = [ — e2/2— leadexp , e2p+leadexp*(l + a/(b+2*n ) ) ] 

27 h = g[0] + t*g[l] 

28 

29 X = var('X') # symbolic Maxima variable 

30 

31 # BEGIN MAIN LOOP: Finds successive coefficients of series h(t) 

32 for i in xrange (2, DEGREE): 

33 h = h + x*t"i + 0(t "(i+2)) 

34 hp = h . derivative () 

35 hpp = hp. derivative () 

36 PEZ = hp*(hpp*t*(l-t))~2 + (hp*(2*h-hp*(2*t-l))+bl*b2*b3*b4)"2 \ 

37 - (hp+bl"2)*(hp+b2"2)*(hp+b3"2)*(hp+b4"2) 

38 LHS = PEZ[i](X) # Substitute x=X in p_i(x) 

39 solucion = solve (LHS==0, X) # Find root of p_i (X) 

40 g . append (QQ( solucion [0 ]. rhs ()) ) # The root becomes 

41 # the next coefficient 

42 # Now reconstruct h up to degree i using the newly— found g[i] 

43 h = g[0] + g[l]*t+g[2]*t"2 

44 for j in xrange (3,i+l): 

45 h += g [ j ] * t " j 

46 # END OF MAIN LOOP 

47 

48 h = h + 0(t "DEGREE) 

49 

50 hi = (h-h[0])/t 

51 F = hi 

52 F — = e2p+leadexp 

sa F /= t-l+0(t "DEGREE) 

54 F = F. integral () 

55 F = F.expQ 

56 F = F. truncate (DEGREE+1) 

57 

58 htilde = hi — e2p 

59 htilde = htilde . truncate (DEGREE-1) 

60 

6i def E(u): 

62 return leadcoef * u'leadexp * F(u) 

63 def Ep(u) : 

64 return ( htilde (u) / (u—1) + leadexp/u/(l — u) ) * E(u) 

65 

66 def pcummul ( t h e t a ) : 

67 return 1-E((1+ cos ( theta ) ) /2) 

68 def nuf theta ) : 
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69 if theta < l.e— 6 or theta > pi — l.e— 6: 

70 return 

71 t = (1+cos (theta))/2 

72 return sqrt ( t*(l — t ) ) * Ep(t) 
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